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ABSTRACT 


Recent work by Peng (‘'Nonlinear Multiaxial Finite Deformation 
Investigation of Solid Propellants**, S.T.T. Peng, AFRDL . TR-85-036) 
has shown that the use of separable symmetric functions of the principal 
stretches can adequately describe the response of certain propellant 
materials and, further, that a data reduction scheme given by him 
gives a convenient way of obtaining the values of the functions from 
experimental data. Based on Peng f s representation of the energy, a 
computational scheme has been developed that allows finite element 
analysis of boundary value problems of arbitrary shape and loading. 

The computational procedure has been implemental in a three-dimen- 
sional finite element code, TEXLESP-S, which is documented in this 


report . 
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1 . 0 INTRODUCTION 

The finite element analysis of large elastic strain 
problems has been reported by many researchers, see for 
example, Oden [1]* Cescotto and Fonder [2], Aly [3], 

Miller [4] Hubbitt et al. [5] and Becker et al. [6]. The 
constitutive relations employed in all of these works is a 
strain energy density which is taken as a function of the 
principal invariants of strain (technically, of the left 
Cauchy-Green deformation tensor ) . No doubt the fondness 
exhibited toward this approach dates from the early work 
of Rivlin, in which it was shown that such forms of the 
energy function are acceptable for any isotropic hyper- 
elastic material. In fact, most published reports of 
hyperelastic solutions employ a very simple form of the 
energy that was introduced by Mooney and has come to be 
known as the Mooney-Rivlin energy function. It is well 
known that this form of the energy describes the behavior 
of real rubber over only a very limitted range of deforma- 
tions, but, due to its simple form (only two constants are 
required to specify the Mooney material) or, perhaps, due 
to the association with Rivlins fundamental work in rubber 
elasticity, the Mooney-Rivlin function continues to be 
used extensively. It is worth noting that the simplicity 
of form, while extremely valuable in the construction of 
analytical solutions to boundary value problems, is of but 
minimal significance in finite element work. 
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Viewed in the large, the stress analysis of any 
problem involves not only the finite element calculations 
used to satisfy equilibrium and boundary conditions but 
also the determination of a constitutive relation appli- 
cable to the material. 

The characterization problem includes determining, 
from experimentally obtained data, both the mathematical 
form of the energy function and the value of the parame- 
ters ("material constants") embedded in these forms. 

While theoretically convenient, the strain invariants 
are far from ideal choices of constitutive variables when 
viewed from the point of view of experimental determina- 
tion of the constitutive form. Literally dozens of forms 
of the energy function, with strain invariants as argu- 
ments, have been proposed in the literature. The ration- 
ales for these are empirical rather than physical - 
material science considerations offer no guide to the 
form. The most common form is a polynomial in the strain 
invariants with the coefficients determined by least 
squares fitting to whatever data are available. Since any 
reasonable function can be approximated by a polynomial, 
so the argument goes, this simplistic approach should, in 
principle, be adequate. Unfortunately, polynomial inter- 
polation invariably produces oscillatory behavior in the 
function and, especially, in its derivatives. Since it is 
the first and second derivatives of the energy that are 
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required in any computational procedure, the high order 
polynomials required to give adequate fits to real material 
behavior over large ranges of deformation do not work well 
in practice. 

The use of principal stretches as arguments of the 
energy function has been proposed by, for example Valanis 
and Landell [7], Ogden [8] and Peng [9]. When the 
stretches are chosen as constitutive variables, an assump- 
tion can be made that reduces the difficulty of the char- 
acterization problem. Using the assumption of separabil- 
ity, Ogden, for example, has been able to fit a wide range 
of rubber deformations using fairly simple mathematical 
forms of the energy function. Ogden's characterization, 
however, requires the choice of some parameters on an 
intuitive basis. 

The approach to material characterization taken here 
differs from previously reported approaches in that no 
mathematical form of the energy function is postulated. 
Since only the values of the derivatives of the energy and 
not form of the function, are needed for computation, we 
avoid the restrictions inherent in the choice of a 
particular mathematical form. Working directly with 
interpolated experimental data we calculate points on a 
curve that defines the material response. These points 
are supplied as input data to the finite element analysis 
code. Our procedure, thereby, provides the most direct 
and straight forward use of experimental data for the 
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solution of boundary value problems. The details of this 
procedure are given in Section 2. 

The finite element code, TEXLESP-S developed in the 
present effort utilizes components of the code TEXGAP3D, 
which is described in reference [10]. The modelling 
capability is identical to that of TEXGAP3D and almost all 
of the data are identical. The hyperelastic calculations 
are based on proceedures developed by Aly in reference 
[3] , but modified to allow the use of the separable energy 
function of principal stretches. The use of TEXLESP-S is 
described in Section 3. 

2.0 THEORETICAL DEVELOPMENT 

The theory of large elastic deformations was essen- 
tially set forth by Rivlin in the late 1940 's. Notations 
have changed and the use of finite element techniques has 
shifted the emphasis but the theoretical foundations are 
unchanged. The major change in emphasis is the use of 
variational principles (principle of virtual work, for 
example) to satisfy equilibrium requirements and the 
accompanying use of Lagrange multiplier methods to accom- 
modate the incompressibility condition. 

In the following sections we review the notation and 
description of deformation, Section 2.1; the constitutive 
relations on which this work is based. Section 2.2; the 
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variational formulation of the boundary value problems. 
Section 2.3 and the finite element implementation of the 
theory. Section 2.4. 


2.1 Kinematics of Deformation 

Let X denote the position of a material particle 
in the undeformed (and unstressed) reference configuration 
of a body whose deformation is to be studied. At some 
later time, t , as the body is deformed by applied loads, 
the particle originally at X moves to a new position 
x . The problem to be solved is the determination for a 
fixed time, the value of x for every particle, X , in 
the body consistent with the requirement that each part of 
the body be in equilibrium with applied loads. We con- 
sider only slow variation in the loads, so that inertia is 
not important, and elastic behavior, so that the history 
of deformation is not important. Thus, time appears in 
all equations only as a parameter identifying which set of 
loads and deformations are being studied. It will be 
convenient to replace time, then, with another parameter, 
p , which we shall call the load factor. The load factor 
will vary, under control of the analyst, from a value of 
zero in the reference configuration through whatever 
positive values are of interest. When the load factor has 
a value of unity, the loads acting on the body will have 
what we shall call their nominal values. We note that 


since all loads are scaled by the same parameter we are 
restricting attention to proportional loading. 
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The motion of the body is described mathematically by 
the function (unknown, a priori ) 

x » x (X, p) (1) 

from which we shall generally omit the dependence on load 
factor p . 

In an elastic material, the stress at a particle X 
depends only on the difference between the shape of an 
infinitesimally small portion of the body around X and 
the shape of this same portion in the reference configura- 
tion. If dS denotes a line element in the reference 
configuration (at particle X) and ds denotes the de- 
formed configuration of this element, then the relation 
between these is given by 

ds = F dS (2) 

where the deformation gradient F is calculated as 
3x 

? - if < 3 > 

The determinant of the deformation gradient, F , 
plays an important role in the kinematics of deformation. 
Direct calculation shows that the ratio of the volume 
contained in a deformed material region dv to the volume 
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of the same material in the reference configuration dV 
is given by 

5? = det M = J (4) 

This determinant, denoted by J , is often callled the 
Jacobian determinant of the deformation or, simply, the 
Jacobian. Clearly, for all physically acceptable defor- 
mations J > 0 and for isochoric (i.e., volume preserv- 
ing) deformations as must occur in incompressible mater- 
ials, J = 1 . 

Since the change in shape of the infinitesimal part 
of the body surrounding the particle X is completely 
determined by the knowledge of the changes of all the line 
segments emanating from the particle, it is clear that F 
contains all of the information required. In fact, F 
contains not only the description of changes in shape of 
the neighborhood of X , but also changes in its orienta- 
tion, since a rigid body rotation of the neighborhood 
would change each dS into a ds with the same length 
but a different direction. A useful theorem of linear 
algebra (the polar decomposition theorem) assures us that 
F can be factored in such a way that the rigid body 
rotation and stretching parts of the deformation are 
separate , Thus , 


F = RU 


(5) 
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Since rigid body rotation does not affect the stress, 
we are interested in the stretch tensor O . Although the 
polar decomposition theorem assures the existence of U , 
it does not offer a convenient way to determine it. But 
since R represents a rigid rotation, R = R , the 
tensor C , defined by (6) , depends only on the stretches. 

C = F T F = U T R T RU = U T U = U 2 (6) 

Clearly, C is a useful measure of the stress producing 
deformation of the neighborhood of the material surround- 
ing the particle X . Another measure of this deformation 
(not used here) is the Green strain tensor, defined as 
E = -|(C - I) 

The tensor C is symmetric, and, therefore, has real 
principal values. If these values are denoted, say, y^ , 
i = 1, 2, 3, then the characteristic equation of C is 

v 3 - IjV 2 + I 2 y 2 - J 2 = 0 (7) 

2 

In (7), Ij, I 2 and J are the three principal invari- 
ants of C . The third of these is written here as 
2 

I^ = J emphasizing the fact that it is equal to the 
square of the Jacobian. 

Since J > 0 , the tensor C is positive definite, 

i.e., y^ > 0 i = 1, 2, 3 . Recalling the last definition 

2 

in (6) , we denote the eigenvalues of C by X rather 
than y . The interpretation of the X A is that they are 
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the principal values of the stretch tensor U or, simply, 
the principal stretches. Rewriting (7) as 

(X 2 ) 3 - I 1 (X 2 ) 2 + I 2 (X 2 ) - I 3 * 0 (8) 

clearly shows that a knowledge of the three principal 
invariants 1^, I 2 , and implies a knowledge of the 

three principal stretches X ^ , X 2 , X^ and vice versa. 
Thus, either set of quantities can, in principle, be used 
to describe the stress-producing aspect of a deformation. 

The geometric interpretation of the tensors F , R 
and U and of the principal stretches are shown, for a 
two-dimensional case, in Figure 1. The circular neigh- 
borhood, N , surrounding X is shown in Figure la. Three 
line elements through X are shown. These are labeled 1, 
2 and 3. Lines 1 and 2 are in the principal directions of 
the stretch tensor U , while 3 is simply another, arbi- 
trarily chosen, line segment. The deformed configuration 
of the neighborhood, n , is depicted in Figure lb. We 
note that, in general, all of the line elements in N 
have been rotated and stretched (or compressed) . This is 
the result of the deformation described by (2) , and we say 
that F carries N into n . If the diameter of N is 
taken as unity, then the lengths of linear segments 1 and 
2 in configuration n (Fig. lb) are X^ and X 2 . These 
are the values of stretch (ratio of deformed to undeformed 
length) for those line elements that get stretched the 
most and least of all line segments through the point X . 
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Figure lc shows the configuration of neighborhood N 
that would be produced purely by the stretching part of 
F , had no rigid body rotation occurred. If we denote 
this n , we say that U carries N into n while R 
carries n into n . We note that there is no rotation 
of the principal line elements, 1 and 2, produced by 0 
and that there are no length changes of any line elements 
produced by R . 

2. 2 Constitutive Equations 

We consider only isotropic, incompressible hyper- 
elastic materials. For this class of materials, the 
stresses are determined, to within a hydrostatic pressure, 
from a scalar function called the strain energy density 
function or, simply, the energy function. In most of the 
published work on rubber elasticity, the energy is written 
as a function of the principal invariants as 

U = Odj, I 2 , I 3 ) (9) 

The Cauchy stress (traction per unit of area in the de- 
formed configuration) is, in terms of the function U(I^, 

X 2' I 3 ) 

a = 2[(U, 1 + IjU^JB - U, 2 B 2 ] + pi (10) 

where U ^ = -|y- 
T 

B = FF 

I is the identity 

p is the hydrostatic pressure 
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In this work, we choose to express the energy func- 
tion, not explicitly as a function of the invariants, but 
rather, following Valanis and Landel, as a separable 
function of the principal stretches. Thus, we write 

U = w(X x ) + wU 2 ) + w(X 3 ) (11) 

The Cauchy stress for an incompressible material for 
which (11) holds is given by 

3 i i 

a = l A. w' (A.) n 1 * n + pi (12) 

i=l 1 1 * 

In (12) , n 1 is the unit vector in the direction of 
the principal stretch A^ . 

For the solution of the equilibrium equations (as 
will be seen in section 2.4), we shall need the first and 
second derivatives of the energy function with respect to 
the invariants. If we were using the form (9) , this cal- 
culation could be made easily and explicitly, but when 
using (11) as our constitutive assumption, we proceed as 
follows . 

Let 

D ’k 5 and “'kH 5 31^17 < 13 > 

denote the derivatives of the energy with respect to the 

2 

three invariants 1^, I 2 and 1^ = J . Also, let 

3 Z k 3 2j k 

J k,i H TT and J k,ij 5 7X73X7 


( 14 ) 
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denote the derivatives of the invariants with respect to 

(14) are 


(15) 


i=j 
(16) 
i=j 

Equating (9). and (11) and differentiating the result 
with respect to A ^ gives three equations for the 
determination of U,^ 

U ,k X k,i = w ' (X i } * - 1* 2, 3 (17) 

We note that for known values of the A^ the coefficients 
I. . in (17) are evaluated using (16) . For a given function 
w (A) we obtain the derivative w' (A^) by the process 
described below. Thus, equations (17) are three linear 
algebraic equations that are solved for the required 
values of U . . 


the principal stretches. The derivatives in 

readily calculated from 

2 2 2 
r l - X 1 + X 2 + 3 

X — A^A^ + A 2 A 2 + A 2 A 2 
I 2 A 1 A 2 + A 2 A 3 + A 1 A 3 


I 3 5 J‘ * X 1 X 2 X 3 


and the results can be written concisely as 


J i,i - 2X i 


I, . . ■ 26 . . 
1/13 13 


I 2 ,i - 2X i (I l " X P 


‘2,ij 


f 2(I i - 


i 4 X i A i 


J ,i ■ J/X i 


'12 


J/X i X i 
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The second derivatives U ^ are obtained in a way 
similar to that just described. Differentiating (17) this 
time with respect to X^ we obtain 


U ,k£ J k,i * 1 , j + °,k *k,ij 

irj = 1, 


(wMX.jOfj 
2 , 3 


(18) 


Equations (18) are six linear algebraic equations for the 
determination of the values of U ^ . We note that the 
values of the first derivatives U ^ must be found by 
solving (17) before equations (18) can be formed. 

When two (or three) of the principal stretches are 
identical, then (17) and (18) are singular, and special 
care must be taken with the degenerate systems. This is 
an exceptional situation that almost never occurs in the 
solution of boundary value problems. TEXLESP-S contains 
adequate provisions to detect and to deal with the 
singular cases. 

In our implementation, values of the function w' (X) 
are input at several values of X . A typical curve 
defined by these data is in Figure 57 of Peng's report 
[9] . Values of w' must be calculated by interpolation 
for use in (17) and values of w" must be calculated for 
use in (18) . The functional values of w' must be smooth 
so that w" will be continuous for all values of X or 
else the tangent stiffness matrix (see eqs. (7) and (10) 
in Section 2.4) cannot be calculated. To assure smooth 
variation of /' with respect to X , we fit the data 
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with cubic splines (using Hermite polynomial interpolating 
functions between data points) . 

2.3 Finite Element Formulation 

The finite element formulation is based on the princi- 
ple of virtual work, which is written as 

J 6U dV * | f • fix dV + | t*6x dS (19) 

ft fl 3 Q 

In (19), the notation is defined as follows: 

6U is the variation of the energy function with 
respect to the current position x 
f is the vector of body force densities 

t is the vector of surface tractions 

$x is the variation of current position 

ft is the interior of the body in the reference 

configuration 

3ft is the boundary of the body in the reference 
configuration. 

We note that the boundary of the body is, in general 
composed of a portion, say 3ft t , on which the tractions, 
t , are given and a portion, 3ft , on which the deformed 
position fo the material is given (and on which the 
traction is unknown but on which the variation of posi- 
tion, $x , vanishes) . Thus, the region of integration of 
the surface integral can be taken simply as 3ft^ . 
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Satisfaction of (19) for all suitably smooth <5x 
guarantees equilibrium (in a variational or weak sense) . 

In order that the deformation be isochoric (i.e. , that the 
incompressibility of the material be satisfied) , we must 
add to (19) the condition that J - 1 = 0 . This is 
accomplished using Lagrange's method of multipliers. The 
result is the modified principle of virtual work 

| 6(u + p (J-l) ) dV = | f • 6xdv + | ffixdS (20) 

Here, p, the Lagrange multiplier, is the hydrostatic 
pressure occurring in the constitutive equation (10) . We 
rewrite the left-hand member (20) as 

| <5u dV (21) 

and note that the variation is with respect to p as well 
as with respect to x and that the function now depends 
on J . 

The body is represented by an assemblage of finite 
elements. The geometry of each element is defined by the 
coordinates of the nodal points connected to it and by a 
set of functions called shape functions, denoted by 
N°(s^) i = 1, 2, 3. The parameters s^ are called the 
parametric coordinates of the material particles contained 
in an element and can be thought of as describing the 
location of a point in a master element. Typically, the 
s^ range over a simple interval such as -1 S & 1 . 
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In TEXLESP-S, the element library consists of 20 node 
bricks, 15 node triangular prisms and 10 node tetrahedra. 
The details of these elements can be found in, for 
example, the book by Zienciewicz [11] or the TEXGAP3D 
reference manual [10]. The shape functions in each type 
of element are quadratic functions of the s^. 

All integrals in the virtual work statement and its 
consequences are evaluated by numerical quadrature 
element -by-element . Thus, all calculations required in 
the following development are performed at certain fixed 
points within an element called integration points. Once 
the finite element discretization is made, there are no 
longer any functions of position (either known or unknown) 
and the only undetermined parameters are the nodal point 
values of position, x , and pressure, p . 

Explicitly, the reference position of a point, the 
deformed position and the pressure at a point are given by 
the following 

X i = N ° X i i = 1/ 2, 3 

x i = N ° x i a * 1, 2, . . . Node (22) 

p = M 8 p 6 8 » 1, 2, 3, 4 

Similarly, the variation of the components of x and of 

p are given by 
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Sx L = N° 6x“ 

6p = M 6 6p e (23) 

For convenience in the following, we will let u stand 

for the entire set of nodal quantities, including both x? 
0 

and p where a ranges over all nodes in the model, and 
8 ranges from 1 to 4 for each element of the model. 
Similarly, other quantities with a sub-tilde, e.g., I , 
will denote vectors with a corresponding number of com- 
ponents. Double sub-tildes will denote square matrices of 
compatible size. 

When the finite element discretizations (22) and (23) 
are substituted into the modified principle of virtual 
work, the results can be written as 

6u T I * 6u T F (24) 

in which 

j - Is | 5av 

* SI 

F = I f *NdV + | t'NdS 

Si dSi 

The equilibrium equation (24) , which is really a vector 
equation with one component for each degree of freedom in 
the problem, is nonlinear since I = I(u, p) is highly 
nonlinear. The solution of this system of equations is 
accomplished by means of incremental loading combined with 


20 . 


Newton iteration. As noted in Section 2.1, we consider 

a 

proportional loading so that we can write f = pf and 
t = pt for the applied loads. Using these conventions in 
(25) , we obtain the obvious modification of (24) 

I = pF (26) 

for the equilibrium equations. The incremental loading 
procedure consists of proceeding from a value of p at 
which (26) is satisfied to an incremented value, say 
p + Ap , and attempting to satisfy (26) again. The first 
step in this consists of a linearization of (26) about the 
current value, i.e., 

31 (u) 

I(u + Au) = I (u) + ■ I" - "- Au 


from which we obtain as the incremental equation 
31 (u) 



Au = 


(p + Ap) F - I(U) 


(27) 


This set of linear equations has as its coefficient matrix 
the tangent stiffness K defined as 



(27) 


We note that it is the evaluation of this matrix and of 
the vector I (u) on the right-hand side of (26) that 
requires the bulk of the calculations in TEXLESP-S. 
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In the first step of solving for a new load increment 
K and I are evaluated at the previous configuration and 

mm +* 

the change in position Au is found by solving ( 26 ) . 

This solution is accomplished in TEXLESP-S using a frontal 
elimination routine. In general, the incremented state 
u 1 = u + Au will not satisfy equilibrium, that is to 
say, the residual 

R 1 = (P + Ap ) F - Ilu 1 ) 

will not be zero. To improve on the solution, we linear- 
ize again — this time about the configuration u 1 . This 
process is 

K 1 du 1+1 = R 1 



+ du 


i+1 


and is repeated until the correction du becomes less 
than some acceptable tolerance. 

In practice, the user of TEXLESP-S specifies a 
sequence of load factors, say p^ , p 2 , P3 ••• P n > and a 
tolerance for the convergence of the du . For each load 
step, the solution, u , is saved on file. Post-process- 
ing, i.e., stress calculation, printing and plotting, can 
be done for any converged load step. The solution proce- 
dure can be continued by restarting the code and specify- 
ing additional load factors. 
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The details of the calculation of I and K are 
contained in the Appendix. 


3.0 TEXLESP-S USERS GUIDE 

The finite element code TEXLESP-S solves equilibrium 
problems for three-dimensional incompressible hyperelastic 
bodies. The constitutive equations are given by specify- 
ing the energy function as either a polynomial in the 
strain invariants or a separable function of the principal 
stretches. In the latter case, the material data are 
input as points on the w' versus x curve. 

The preprocessing and post-processing functions of 
TEXLESP-S have been adapted from the code TEXGAP3D for the 
solution of three-dimensional linear elastic problems. 
Insofar as possible, the differences between these codes 
have been made transparent to the user. That is to say, 
all modeling and post-processing data are identical in the 
two codes. Only the material property descriptions and 
commands directing the solution procedures are different. 
Consequently, the TEXGAP3D users manual can be used, ver- 
batim, for all descriptions of input formats, mesh genera- 
tion commands, element definitions, boundary condition 
specifications and post-processing commands. 

This users' guide contains an overview of the struc- 
ture of the data deck (Section 3.1); a description of 
those data commands that are different from TEXGAP3D data 
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(section 3.2); and a set of example problems which 
illustrate the use of the code (Section 3.3). 

3. 1 Structure of Code and Data Deck 

Data deck structures for TEXLESP-S and TEXGAP3D are 
essentially the same. In each code, there are data 
required to 

a) set up the model 

b) solve the problem 

c) post-process the solution. 

While these are typically all contained in a TEXGAP3D run, 
they may occur in separate runs for TEXLESP-S. Using the 
restart capability in TEXLESP-S is often desirable and 
requires some understanding of the code structure. 

Figure 2 shows a large-scale flow diagram of 
TEXLESP-S. It can be noted that any run begins with the 
main routine. From the main routine, control can be 
directed to 

a) SETUP - to define a model 

b) SOLVE - to calculate solutions for various 
load increments 

c) POST - to calculate print and plot 
stresses, strains and/or displacements 


d) RESTART - to resume execution of problems 
begun on a previous run. 
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It will be noted that the file TAPE12 is used by all 
modules. This file contains all data that describe the 
model, i.e., nodal point coordinates, element definitions 
and boundary condition specifications, as well as the 
solutions calculated at each load step. Data are written 
to TAPE12 by SETUP and by SOLVE. Any run which calls 
these modules modifies TAPE12 and, if subsequent process- 
ing of the job is to be performed, the updated version of 
this file must be saved. 

TAPE 18 is a file that contains the material data for 
hyperelastic materials whose energy function is written as 
a separable function of the principal stretches. These 
data are read by TEXLESP-S in the SOLVE and POST modules 
but are not modified. This file must be attached to the 
job when SOLVE or POST are to be used. The format of 
TAPE18 is such that the following statements can be used 
to read it. 

READ (18 , 2010)N 

READ (18,2020) (ALAM (I) , I = 1,N) 

READ (18,2020) (WP(I) , I = 1,N) 

2010 FORMAT (15) 

2020 FORMAT (6F1 1.3) 

The data read by these statements are: 

N the number of points on the w' vs X curve 
ALAM (I) values of X , in ascending order 


WP (I) values of w' 
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Data Deck Structure . The structure of data decks for 
TEXLESP-S will vary according to the nature of the job 
being processed. In Figure 3, we show the most general 
structure, i.e., that which is used when a problem is to 
be solved in a single run. The data cards that are given 
explicitly are similar for all such jobs. Data that are 
problem-dependent are not given explicitly in Figure 3. 

Figure 4 shows the structure of a data deck that 
would be used to continue a previously started solution 
and to perform post-processing on it. 

3. 2 Description of TEXLESP-S Commands 

Only those data cards that differ from the corre- 
sponding TEXGAP3D data are described here. For complete 
descriptions of the other data, see the TEXGAP3D Users' 
Manual [10]. The following descriptions explain all data 
cards that differ from those of TEXLESP3D. 

3.2.1 Integration control card . 

SETINT 

This card, which should occur before the SETUP card 
changes the default number of integration points from 14 
to 8. A significant saving in computer time can be 
realized by the use of the SETINT card. 
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3.2.2 Material property definition cards 


RUBBERS , name , no . 


This card defines the material, "name," as material 
number, "no" and specifies that it is a rubber elastic 
material whose energy function is defined by the data on 
TAPE 18 . 

Note that only one such material can be used in a 
problem. 


RUBBER I, name, no, n, CIO, C01, C20, Cll, C02 

C30 , C21 , C12, C03 


This card defines material, name, as material number 
"no." The energy function for the material is defined in 
terms of the strain invariants in the following way, 
depending on the value of n . 

n = 1 V 1 = C10(I 1 -3) + C01(I 2 -3) 

n = 2 U 2 = u i + C20 (Ij-3) 2 + ClMIj-3) d 2 “3) + C20(I 2 »3) 2 

n = 3 U 3 * U 2 + C30(I 1 -3) 3 + C21 (I 1 -3) 2 (I 2 -3) 

+ C21(I 1 -3) (I 2 -3) 2 + (C03(I 2 -3) 3 

These are the only two material types available in 


TEXLESP-S 
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3.2.3 Grid Generation 

The grid generation commands in TEXLESP-S are 
identical to those in TEXGAP3D. See Section 2.4 of the 
TEXGAP3D Users' Manual for descriptions of these commands. 

3.2.4 Element Definition and Boundary Condition 
Specification 

The element library and element generation commands 
are identical in TEXLESP-S and TEXGAP3D except that the 
singular WEDGE elements are not used in TEXLESP-S. 

All of the boundary condition commands in TEXGAP3D 
are available in TEXLESP-S. In addition, a CLAMP command, 
described below, has been added to TEXLESP-S. 

The element and boundary condition specification 
(with the exceptions noted above) are described in Section 
2.5 of the TEXGAP3D users' manual. 


BC, CLAMP, i, j, k, nside, nset, vx, vy, vz 


This boundary condition specifies that all of the 
nodes on side, "nside" of element "i, j, k" have specified 
displacements. The values of these displacements are vx, 
vy and vz in the x, y and z directions 
respectively. Each different set of vx, vy, vz is given 
a number, "nset," the first time it is used. When a given 
set of applied displacements is used subsequently only the 
set number (n set) need be repeated. See example data set 
1 for use of the CLAMP command. 
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3.2.5 Solution commands 

In TEXGAP3D, the command, SOLVE, has no parameters, 
since only a linear solution is produced. In TEXLESP-S, 
there are parameters on the SOLVE card that control the 
incremental loading and Newton iteration processes. 


SOLVE, jprint, iter, tol, p 2 , p 3 , ..., p m 


The SOLVE command causes the solution of the 
equilibrium equations to be performed using the sequence 
of load factors p^, p 2 , p^, . .., P m . Each solution is 
recorded on TAPE12. If no previous solutions have been 
stored on TAPE12, then these are the first n solutions. 
On the other hand, if a previous run of the problem has 
produced solutions and the current run is a RESTART run, 
then the sequence of solutions produced is recorded fol- 
lowing the solution at which the restart was made. See 
example problem 2 for an illustration of this. 

The parameter, jprint, controls the amount of print 
produced during the solution. Increasing values of jprint 
produce increasing amounts of output. A value of 1 pro- 
duces the convergence summary table and is recommended. 

The maximum number of iterations allowed in a load 
step is equal to "iter." The default value, which is 10, 
is recommended. 


I 
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The parameter, "tol," defines the tolerance to which 
the changes in displacement are compared when determining 

whether the Newton iterations have converged. The default 

-4 

value is 10 which is small enough for nearly all 
purposes. The convergence check determines that 

du^ < tol * |UjJ 

for each component of displacement in the model. Here, 
u^ is the value of the displacement component and du^ 
is the correction to that displacement produced by the 
Newton iteration. 

A summary of the convergence tests is printed for 
each load step. If the iterations should fail to converge 
at any load step, all of the solutions up through the 
previous load step are stored on TAPE 12 so that restarting 
can be done from that point. 

3.2.6 Post Processing 

The post processing in TEXLESP-S uses commands that, 
with two exceptions, are the same as those used in 
TEXGAP3D. All plotting and stress calculation commands 
are as described in Section 2.7 of the TEXGAP3D users' 
manual 


POST, nstep 
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The POST command in TEXLESP-S contains the parameter 
"nstep" which specifies the load step number of the 
solution to be processed. 


REACTION, xmin, ymin, zmin, xmax, ymax, zmax, 
imin, jmin, kmin, imax, jmax, kmax 


The REACTION command causes the reactions at all 
nodes within the specified (x,y,z) bounds and the speci- 
fied (i,j,k) bounds to be printed. Either or both sets of 
bounds can be omitted. If no bounds are specified, then 
all nodes in the model are processed. 

The reactions that are printed are actually the un- 
balanced nodal point forces (i.e., the difference between 
specified applied loads and the calculated internal 
forces). When equilibrium is exactly satisfied at a node, 
the reactions will be equal to zero, unless the node has 
had displacements specified. When displacements are 
specified, then the actual reaction furnished by the 
support on the model is the calculated value. 

The vector sum of all of the reactions in the 
specified part of the model is also printed. The REACTION 
command is, thus, very useful in determining load-deflec- 
tion data. 
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3.3 Example Problems 

The example problems presented in this section are 
intended to illustrate the application of TEXLESP-S in a 
typical application. Descriptions of the input data for 
several runs associated with the problem are described 
below. 

The problem to be analyzed is the stretching of a 
thin biaxial test specimen with a rail cross section. The 
specimen is shown in Figure 5. This is the same specimen 
for which results are given in Peng's report [9] (see Fig. 
A- 2 in [9]) . 

The separate runs that are described below are: 

Run 1 - Generation and plotting of grid 
Run 2 - Analysis of stretching up to 20% 

Run 3 - Analysis of stretching from 20% to 90% 
Run 4 - Calculation of reactions 
Run 5 - Stress calculation and deformed shape 
plotting. 

It is desirable, although not necessary, to perform 
nonlinear analysis in several stages using the restart 
capability of TEXLESP-S. The use of several separate runs 
allows the user to check intermediate results to make 
sure, for example, that the grid is error free and that 
convergence of the nonlinear solution is taking place. 

For runs 2 through 6, TAPE 12 from the prior runs and 
TAPE 18 containing material properties are required. 
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Run 1: Generation and Plotting of a Grid 

The first step in generating a grid for TEXLESP-S is 
to sketch the grid and to assign (I, J, K) node numbers. 
Figure 6 shows such a sketch of one-eighth of the speci- 
men. In this grid, the I, J, K directions are in the x, 

y, 2 directions respectively. Since the elements are 
ordered, for solution purposes, by varying I then J 

then K it is especially important to let I vary in the 

direction of fewest elements (one or two in this case) and 
K in the "long" direction (5 elements in this case). It 
will be noted that several nodes have been identified on 
the sketch. 

The data deck for this run is given in Figure 7. 
Following the title card, there is a SETINT command 
(Section 3.2.1). Then the SETUP command calls the 
pre-processing overlay. 

The material model in this problem is given by the 
w' versus A data which are contained on TAPE18. The 
RUBBERS card specifies this and identifies the material as 
material number 1. 

The grid generation proceeds by defining all of the 
nodes on the plane z = 0 which is also K = 1 . FACE, 
ARC, POINT and CONNECT commands are used to define these 
points. After all of the K = 1 points have been 
defined, the NORMAL command is used to define other K 
planes. After all nodes have been defined, the element 
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definition begins. The elements in this problem are 
BRICKH and PRISMH elements. The looping features of 
TEXLESP-S are used in defining the elements. 

The boundary conditions on the one-eighth model are 

(a) Symmetry (SLOPE) on the plane y = 0 which is face 5 
of elements (1, 1, K) . 

(b) Symmetry (SLOPE) on the plane x = 0 which is face 4 
of elements (1, 1, K) and (1, 3, K) . 

(c) Symmetry (SLOPE) on the plane z = 3.75 which is face 
3 of elements (1, 1, 9) and (1, 3, 9) and face 4 of 
the PRISMH, (3, 3, 9). 

(d) All displacements specified (CLAMP) on the plane 

y = .75 . This is face 2 of elements (1, 3, K) and 
face 1 of the PRISMH elements (3, 3, K) . The values 
of displacements specified in the CLAMP command 
(Section 3.2.4) are u = 0, v = .75, w = 0. Note 
that this corresponds to a stretch of 100%. This 
value would be achieved by a load factor of 1.0. 

After the elements and boundary conditions have been 
specified, the preprocessing phase is ended. 

In order to examine the grid generation results 
before attempting an analysis, this run ends with a call 
to the post processor to plot the grid. The plot produced 
by this call is shown in Figure 8. In addition to the 
plot, a large amount of printed output is generated, from 
which detailed checking of nodal coordinates, element 
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connectivities and boundary condition specification can be 
done. 

The data generated by this run are written on TAPE 12 
so that subsequent runs need not repeat the grid genera- 
tion process. 

Run 2 ; Analysis of Stretching Up to 20% 

Although final results may be desired only for a 
couple of values of stretch, we choose to increment the 
load rather slowly at first. This strategy has the 
benefit of producing more complete load-deflection data. 
More significantly, it is often the case that smaller load 
steps are necessary in the early stages of loading in 
order to achieve convergence. In the present example, for 
instance, a load increment from 0 to 20% will not con- 
verge while steps from 15% to 20% as well as from 20% to 
40% do converge. 

Examination of the load-deflection curve for this 
example. Figure 9, indicates a rapid change in behavior 
around a stretch of 10%. It is this behavior that is 
reflected in the need for initially small load steps. 

The data set for this run is shown in Figure 10. 

This run requires attachment of TAPE 12 from Run 1 and 
TAPE 18 containing the material data. The only output from 
this run describes the convergence of each load step. The 
resulting solutions are accumulated on TAPE12. 
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Run 3: Analysis of Stretching from 20% to 80% 

After having determined, from the results of Run 2, 
that the solution is converging satisfactorily, this run 
is used to continue the solution for subsequent load 
steps. Note that the restart card contains the specifica- 
tion of the load step beyond which the next SOLVE command 
will calculate solutions. TAPE 1 2 and TAPE 18 are required. 

Run 4t Calculation of Reactions 

In this run, the code is restarted and the solutions 
at load steps 1, 2 and 3 are post processed. The REACTION 
command produces node-by-node reactions for all nodes 
within the xyz bounds specified. These are not as useful 
as the total resultants that are the sum of these nodal 
reactions. The region specified in this run includes the 
plane y = .75 but not the plane y *= 0 . Thus, the y 
reaction printed is the total applied load on one fourth 
of the top grip. The total applied load is, of course, 
four times this value. 

Figure 9 shows the load deflection data plotted by a 
separate plotting program. 

Run 5: Stress Calculation and Deformed Shape 

Plotting 

To obtain values of displacement, stress and/or 
strain at various locations in the model, the post 
processing commands that are common to TEXLESP-S and 
TEXGAP3D are used. In this example, we calculate dis- 


36 . 


placements, stresses and strains for load step 3. The 
data set shown in Figure 10 also produces a deflected 
shape plot. 
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Figure 2 . Schematic of TEXLESP-S 
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$TITLE 

SETINT 

SETUP 

(Material property data cards) 

END 

(Nodal point generation data cards) 

END 

(Element definition and boundary condition data 
cards) 

END 

SOLVE 

POST 

(Post processing commands) 

END 

STOP 


Figure 3 . Data Deck Structure for Complete Pun 



$TITLE 

RESTART 

SOLVE 

POST 

(Post processing commands) 

END 

STOP 


Figure 4 . Data Deck Structure 


for Post Processing Run 



W.50 








44 


$RUN 1 GRID GENERATION hOR THIN UNCRACKED BIAXIA 1 T-NS ri - C PC T 

SET I NT 

SETUP, 1 

RUBBERS, 1, 1 

END, MAT 

FACS-L, 1, 1, 1, 1, 3, 3, 1 

0, 0, 0, . 06£5, 0, 0, . 0625, . 312, 0, 0, . 312, 0 
FACE-L, 1, 1, 3, 1, 3, 5, 1 

0, . 312, 0, . 0625, . 312, 0, . 0625, . 75, 0, 0, . 75, 0 

ARC-L, 1,3, 5, 1,5, 5, 1, , . 0625, .75, 0, . 5, . 75, 0 

POINT, 1,4,4, 1,. 1903, .6217,0 

POINT, 1, 4, 1, 1, . 1903, 0, 0 

POINT, 1, 5, 1, 1, . 5, 0, 0 

CONNECT, 1, 4, 1, 1,4, 4, 1 

CONNECT, 1,5, 1, 1,5, 5, 1 

NORMAL, 2, 1, 1, 1,5, 5, 1,. 25 

NORMAL, 6, 1 , 1 , 3, 5, 5, 3, 2. 75 

NORMAL, 2, 1 , 1 , 9, 5, 5, 9, . 75 

END, GRID 

KLOOP, 5, 2 

JLOOP, 2, 2 

BRICKH, 1, 1, 1, 1 

BC, SLOPE, 1, 1, 1,4 

JEND 

PRISMH, 1, 3, 3, 1, 5, 5, 1, 3, 5, 1 

BC, SLOPE, 1, 1, 1, 5 

BC, CLAMP, 1,3, 1,2, 1,0, .75,0 

BC, CLAMP, 3, 3, 1, 1, 1 

KEND 

JLOOP, 2, 2 

BC, SLOPE, 1, 1, 9, 3 

JEND 

BC, SLOPE, 3, 3, 9,4 
END, ELEMENTS 
POST, SETUP 

PLOT, ELEMENTS, , , -1 , -1 , -1 , 2, 2, 5 

ANGLES, 0,0,0 

END 

STOP 


Figure 7. Data Deck for Grid Generation Run 
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Figure 9 . Load Deflection Curve 







4 RUN £ ANALYSIS OF STRETCHING UP TO £0* 

RESTART, SETUP 

SOLVE, 1, , , .05, . 1, . 15, .2 

STOP 


4RUN 3 ANALYSIS OF STRETCHING FROM £0* TO 80% 
RESTART, LAST 
SOLVE, 1, , , . 3, . 6, .8 
STOP 


4 RUN 4 CALCULATION OF REACTIONS 

RESTART 

POST, 1 

REACTION, -1, .7, -1,5, 5, 5 
POST, £ 

REACTION, -1, . 7, . 1, 5, 5, 5 
POST, 3 

REACTION, -1, . 7, -1,5, 5, 5 

END 

STOP 


4 RUN 5 STRESS CALCULATION AND PLOTTING DEFORMED SHAPE 

RESTART 

POST, 3 

BLOCK,,, ,-1,-1, -1,2, 5, 5 
OPTION, £ 

PLOT, DEFORMED, ,-1,-1, -1, 1,3,4, S 

ANGLES, 0,0,0 

END 

STOP 


Figure 10 . Data Decks for Restart Runs 
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APPENDIX 

CALCULATION OF ELEMENT MATRICES 


For a Rivlin polynomial material, the contribution of 
an element to the internal energy of the body is 


I 


? 

u(I.,I 2 ,J) + P(J-l) dv 

n e 


where the energy has been augmented to include the incom- 
pressibility constraint (J-l * 0) . The contribution of 
the element to the right-hand side (or load vector) is 
-61 . The variation in the energy is obtained by varying 
all the nodal displacement variables in the element (x?) 
and the pressure variables (p Y ) . Therefore, 


-61 




+ 


P 




+ M Y (J-l) 6p Y }dv 


where 

A « • 

u, k = yj— (note: I^ denotes J = det(F)) 

k 

M y = finite element interpolators for the pressure. 


The contribution of the element to the tangent stiffness 
matrix is d(6I) and is given by 


361 ,6 , 361 . n 
— fi dx i + — n d P 
3Xj 1 3p n 


d (6 1) 
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31. 31 3 2 I 

{ ^ u 'k£ a ¥ + u 'k ci 6 + P 
fi e K£ 3x“ 3Xj K 3x?3x! 3xpx^ 


f 3I k 

■ I 




+ M n 5x?dp n + M y ■—■ 5p Y dx^}dv . 

3 x® 1 3x7 3 


The first and second derivatives of the energy density 
function (u,^ and u, k£ ) can be obtained explicitly from 
the polynomial form of the energy function for Rivilin 
materials or as outlined in section 2.2 for materials 
defined by a m 1 vs. X curve. 

The invariants are complicated functions of the nodal 

O 

point displacements (x?) and derivatives of the invariants 
are best obtained by considering the derivatives of the 
deformation gradient (F) , the left deformation tensor (B) 
and using the chain rule. Substituting the finite element 
interpolation of the nodal displacements into the defini- 
tion of the deformation gradient gives 


mn 


3x 

ax 


3NJ 

ax. 


n 


m 


N y x y 
n m 


(where N Y = jf- ) 
n 


The first and second derivatives of F are: 
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3F 


3x : 


— = N°$ . 
a n im 


3 2 F. 


mn 


3x?3x® 


= 0 


where 6. is the Kronecker delta function. 
im 

The left deformation tensor is defined by 

B = F F 
mn mp np 

The first and second derivatives of B are 


3B 


mn 


3F 




3F 

F + F 

mp 3x « np 3x a 


mp 


= (F 6 + F 6 )N 

mp m np im p 


3 2 B. 


mn 


3x?3Xj 


(6.6. +6.6. )N a N 6 

' jm in ;jn im p p 


(6.6. + 6 . 6. ) 6 N°N e 

jm in jn im pq p q 


The derivatives of the invariants (I^Ij^J) can 
be determined: 


I, = 6 B 
1 mn mn 


Ih 


3B. 


= 6 


mn 


11111 3x“ 


= 2F.lT 

ip p 


now 
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= N 0 ^ 1 . 

P Pi 


3 2 I. 


3 2 b 


3x?3x? 


= 6 


mn 


11111 3x?3x|j 


= 26 . .6 N°N 6 
13 pq p q 

= hJ. n“n? . 
ijpq P q 


The second invariant is more complicated. 


Let (B) = B B 

mn mr rn 

then 


I 1 = tr (B ) 


I, * 6 B B 
1 mn mr rn 


= B B 
mr rm 


The derivatives of 1^ are : 
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3x® 


3B, 


= 2B_ 


mr 


mr 3x? 


= 2B (F 6 . + F 6 . )N° 

mr mp lr rp im' p 


= 4B . F N 
lm mp p 


3x?3Xj 


3B. 

4 [F — + B 


mp 3x! ' ~ iin avS "P 


3x 


rcP J 


= 4 [F (F. 6 . + F 6 . . ) + B, i. { ]nV 

mp }m mq im P3 P <3 




* 4 [F . F, + 6 . .C + B . . 6 ] N“N 

dp iq ij pq ij pq p q 


where C = F F is the right deformation tensor. 


The second invariant is defined by 

I 2 = I (1 1 " 


and its derivatives are: 
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i 3Ii 3 1. 

1( 2 I 1 i] 

21 1 3x“ 3 xf 


I. (2F. N“) - 2B. F N“ 

1 ip p lm mp p 


2(F. I. - B .F ) N“ 
ip l im mp p 


= 2F (6.1. - B. )N“ 
mp im l im p 


N P B Pi 


» 2 i. 


3x i 3x j 


lh!h + I 

3x|? 3x° 1 3x?3Xj 


1 


2 3x?3Xj 


= [4F . F. + 21,6. .6 - 2(F. F. 

]q ip i id pq dp iq 

+ 6 . .C + B . . 6 ) ]N a N 6 

id pq id pq J p q 


2 [ 2F . F. - F. F. + (1,6. . - B..)6 

ip jq ]p iq '1 13 ID Pq 

- C 6 . . ] N a N 8 

pq i] p q 


= H?. N a N 6 

iDpq P q 


The determinate of F is defined as 
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where is the cofactor of . 

The first derivative of J is 


— - 

3x“ *P ? 


- «; j <rV 

- n“b|. 
p pi 


where F ^ is the inverse of the deformation gradient, 
In order to take the second derivatives of J , it is 
necessary to write the first derivatives explicitly. 


- N ?« p 22 F 33 - P 23 P 32> + 

N 2 (F 23 F 31 “ P 21 F 33 ) + N 3 (F 21 F 32 “ F 22 F 31 ) 


* X?IP 13 P 32 * P 12 P 33> + 

N 2 ,P 11 P 33 - P 13 P 31> + N S< P 12 F 31 ‘ P U P 32 ) 


^ - N ! ,P 12 P 23 * P 13 P 22> + 

N 2 (P 13 P 21 - P 11 P 23> + N 3 ,P U P 22 ‘ P 12 P 21> 


The second derivatives are: 


3 2 J 


3x“aXj 


= 0 


if 
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otherwise , 


3 2 J 

ox°3x 2 


,8 


6, 


,6 


= N^(F 33 N^ - F, 0 N*) + N“ (F 01 N, - 


323 2 313 

+ N° (F N® 
w 3 l 32 W 1 


f 33 nJ) 


P 31 N I> 


3 2 J 


3x°3x 3 


r B 


B, 


,8 


= N“(F 22 N^ - F„N”) + N“ (F 00 N, - 


23 2 


23 1 

+ N a (F N® - 
W 3 '* 21 W 2 


F 21 N I> 


F N®) 
* 22l* 


3 2 J 


3x 2 3x 3 


= N“(F 13 N^ - F 12 N^) + N“(F n N5 - F 13 nJ) 

+ N 3 ,F 12 N ? - P 11 N I> 


or 


3 2 J 


3 2 J 


3x?3x? 


3x?3x| 


i ¥ j 


3 2 J 


3x?3x? 


tt 3 KT a V^ 

H ijpq Vq 


Let R? denote the component of the elemental right- 

th 

hand side vector in the i direction for node a and 
R Y denote the component of the vector corresponding to the 
pressure variable. Then, using the expressions for 
the derivatives of the invariants 


* 1 - - \ j u ’i 5 pi + u >2 § pi + <“'j + p> 5 pi>»; 
n e 

R Y = - j M Y (J - 1) dv 

V 
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Similarly, let K?? denote the component of the tangent 
stiffness matrix coupling the degree of freedom in 
direction i at node o with the degree-of-freedom in 
direction j at node 6 , K? n denote the component 
coupling the n t ^ 1 pressure degree of freedom with the 
degree-of-freedom in direction i at node a , and K Yn 
denote the component coupling the and n^ pressure 

degrees of freedom. The stiffness matrix is then 


K ij 


K° n 


-l 

i 

'I 


[u,. . B k . B £ . + u,, H k . + pH?. ]N a N B dv 
1 'kl pi qj 'k ljpq F ijpq j p q 


M 11 B 3 . N° dv 
» Pi P 


K Yn = 


0 


